	
*****************************************************************		
* set up dataset on which to collect counterfactual holdings for each library 
	
preserve 
		clear 
		set obs 1 
		gen var1=. 
		
		save "$datapath/intermediate/collect_pdouble_income.dta", replace		
restore 
*****************************************************************		
	
	
*****************************************************************		
** Define the program that'll maximize librarian utility

capture program drop find 

program define find
		tempvar Dgp 
		tempvar Dge 
		tempvar sDg 
		tempvar Qp
		tempvar Qe
		tempvar b
		tempvar c
		tempvar d
		tempvar e
		tempvar f
		tempvar ff
		tempvar g
		tempvar con
		tempvar ne
		tempvar np
		tempvar n2

     if "`3'"!=""{
        loc sse = "`3'"
        qui cap drop `sse'
        }
     else {
        tempvar sse
        }
		

			gen double  `Dgp' = (abs(`1'[1,1])*exp($dpvar/(1-$sbvar)))^((1-$sbvar)/(1-$stvar)) 
			gen double `n2' = ($bvar - abs(`1'[1,1])*$ppvar)/$pevar 
			gen double `Dge'  = (abs(`n2')*exp($devar/(1-$sbvar)))^((1-$sbvar)/(1-$stvar))
			gen double `sDg'  = `Dgp' + `Dge' 
			gen double `Qp'  = $mvar*`Dgp'*(`sDg'^(-$stvar))/(1+`sDg'^(1-$stvar))
			gen double `Qe'  = $mvar*`Dge'*(`sDg'^(-$stvar))/(1+`sDg'^(1-$stvar))
			gen double `np' = abs(`1'[1,1])
			gen double `ne' = abs(`n2')
			
			gen double `f' = ($tvarp*`Qp' + $tvare*`Qe')  - 1000000*abs( $bvar - $pevar*abs(`n2') - $ppvar*abs(`1'[1,1]) )
			gen double `con' = $bvar - $pevar*abs(`n2') - $ppvar*abs(`1'[1,1])
			su `f' `con'
			su `np' `ne'
			
			egen double  `sse' = sum(`f')
		
			sca `2' = `sse'
		
		
	end
*****************************************************************	


use $datapath/clean/main_data.dta, clear 

merge m:1 fipsplac fipsst using "$datapath/clean/income_place_data.dta" 
	drop if _merge==2
	drop _merge 

	xtile ino = medinc, n(5)
	drop if medinc==.	

	ren prmatexp Ep
	ren elmatexp Ee
	
	
	forvalues i = 1(1)5{
		foreach f in p e {
			su N`f' if ino==`i'
				local N`f'`i' = r(mean)
			su E`f' if ino==`i'
				local E`f'`i' = r(mean)
				
			local P`f'`i' = `E`f'`i'' / `N`f'`i''
		}
	}
    
	di `Pp1'

*****************************************************************		
** Create the sigma and delta variables from demand estimation

	reshape long Q N, i(popu_lsa year lno ) j(type) s
	
	gen double M = 50*popu_lsa
	gen double s = Q/(N*M) 
	
	drop if s==.
	keep if s>0 & s<1 
	
	gen double de = type=="e"
	egen double QQ = sum(Q), by(lno year)
	
	gen double ls = ln(s) - ln((M - QQ)/M)
 
	gen double lstop = ln(Q/QQ)
	egen double NN = sum(N), by(lno year)

	gen double lstop_inst = ln(N/NN)
	gen double lsbottom = ln(1/N)
	gen double lsoverall = ln(1/NN)
	

	ivreghdfe ls  (lstop =   lstop_inst ) lsbottom   i.year de, absorb(lno) robust

	gen double bstop = _b[lstop]
	gen double bsbottom =_b[lsbottom]
	
	
	gen double del = ls - bstop*lstop - bsbottom*lsbottom   
*****************************************************************		
	

keep if year==2018 
	
	
*****************************************************************		
** prepare variables (derivatives) needed to calculate the thetas
	
	gen double Dg = (N*exp(del/(1-bsbottom)))^((1-bsbottom)/(1-bstop))
	egen double sDg=sum(Dg), by(lno year)
	gen double shat = (1/N)*Dg*(sDg^(-bstop))/(1+sDg^(1-bstop))
	
	gen double N1p = N + .1*(1-de)
	gen double N1e = N + .1*de 
	
	gen double Dg_1p = (N1p*exp(del/(1-bsbottom)))^((1-bsbottom)/(1-bstop))
	egen double sDg_1p=sum(Dg_1p), by(lno year)
	gen double shat_1p = (1/N1p)*Dg_1p*(sDg_1p^(-bstop))/(1+sDg_1p^(1-bstop))

	gen double Dg_1e = (N1e*exp(del/(1-bsbottom)))^((1-bsbottom)/(1-bstop))
	egen double sDg_1e=sum(Dg_1e), by(lno year)
	gen double shat_1e = (1/N1e)*Dg_1e*(sDg_1e^(-bstop))/(1+sDg_1e^(1-bstop))

	
	* note: dqp is the derivs of the two q's wrt p 
		gen double dqp = M*(N1p*shat_1p - N*shat)/.1 
		gen double dqe = M*(N1e*shat_1e - N*shat)/.1 

	
	** Set prices 
		** materials costs calculated as mean expenditure/mean holdings by income quintile

		gen pricep = `Pp1' if ino==1
		replace pricep = `Pp2' if ino==2
		replace pricep = `Pp3' if ino==3
		replace pricep = `Pp4' if ino==4
		replace pricep = `Pp5' if ino==5
		
		gen pricee = `Pe1' if ino==1
		replace pricee = `Pe2' if ino==2
		replace pricee = `Pe3' if ino==3
		replace pricee = `Pe4' if ino==4
		replace pricee = `Pe5' if ino==5

		gen pp = pricep
		gen pe = pricee
*****************************************************************		
		

	
*****************************************************************		
** Calculate thetas based on sigmas (and relative prices)

	preserve 
	
		keep lno type dqp dqe pp pe 
		reshape wide dqp dqe, i(lno ) j(type) string
		rename  dqpp dqp_dnp  
		rename  dqee dqe_dne  
		rename  dqep dqp_dne  
		rename  dqpe dqe_dnp  
		
		gen double det = dqp_dnp*dqe_dne - dqp_dne*dqe_dnp 
		
		gen double thetap = (dqe_dne/det)*pp - (dqe_dnp/det)*pe 
		gen double thetae = -(dqp_dne/det)*pp + (dqp_dnp/det)*pe 
		
		keep lno theta*
		reshape long theta, i(lno) j(type) string 
		drop if theta==. 
		egen double nn=count(theta), by(lno)
		keep if nn==2 
		save theta.dta , replace 
		
	restore 

	merge 1:1 lno type using theta.dta  
	keep if _merge==3
*****************************************************************		
	
	


*****************************************************************		
** Turn it into a wide dataset 	
	gen price = pp
	replace price= pe if type=="e"
	
	keep N theta bstop bsbottom price M del lno type Q ino 
	
	
	reshape wide N theta  price del Q, i(lno bstop bsbottom M) j(type) string 
	
	gen ratio = thetae/thetap 
	
	su ratio, de 
	
	gen double budget = Ne*pricee + Np*pricep
	
*****************************************************************		
	
	
	


*****************************************************************		
*****************************************************************		
** Find optimal mNp and mNe for each library
*****************************************************************		
*****************************************************************		


* randomize the order 
	gsort lno 
	set seed 86
	drawnorm rando 
	gsort -rando 

local mmm =_N  


forvalues j= 1 (1) `mmm' {	


************************************		
** baseline and double Pe


	foreach f in 10 20 { 
	
	preserve 

	    replace pricee =pricee*(`f'/10)
	
		keep if _n==`j'
	
		gl nlist = "beta:A"	
		gl dpvar = "delp"
		gl devar = "dele"  
		gl bvar = "budget"  
		gl stvar = "bstop"  
		gl sbvar = "bsbottom"
		gl ppvar = "pricep"  
		gl pevar = "pricee"  
		gl tvarp = "thetap"  
		gl tvare = "thetae"  
		gl mvar = "M"  
		gl np0var = "Np"  

		su Np 
		local p=r(mean)
		
		su Ne
		local e=r(mean)

				 
		 mat b0 = ( `p' )
		 mat colnames b0 = $nlist
		 
		 amoeba find b0 yout bols . . 1E-10 
		 
*****************************************************************		

		 
		 
*****************************************************************		
** Calculate holdings, quantities and CS based on optimal holdings

		 
		gen mNp = bols[1,1]
		gen mNe = (budget - pricep*mNp)/pricee 
		
		su Np Ne mNp mNe 
		gen factor = `f'
			gen double  Dgp = (mNp*exp(delp/(1-bsbottom)))^((1-bsbottom)/(1-bstop)) 
			gen double Dge  = (mNe*exp(dele/(1-bsbottom)))^((1-bsbottom)/(1-bstop))
			gen double sDg = Dgp + Dge 
			gen double mQp  = M*Dgp*(sDg^(-bstop))/(1+sDg^(1-bstop))
			gen double mQe  = M*Dge*(sDg^(-bstop))/(1+sDg^(1-bstop))
			gen double CS  = M*log(1+sDg^(1-bstop))
		
*****************************************************************		
	

*****************************************************************		
** save the variables we care about

		keep Np Ne mNp mNe lno factor Qe Qp  mQe mQp CS lno M dele delp thetae thetap ino

		gen raise ="Pe"

		sleep 500
		append using "$datapath/intermediate/collect_pdouble_income.dta" 
		sleep 1000
		capture save   "$datapath/intermediate/collect_pdouble_income.dta", replace  

	restore 
		
	}

**************************************	
	
	

************************************		
** double Pp
	
	preserve 

		replace pricep = 2*pricep
	
		keep if _n==`j'
	
		gl nlist = "beta:A"	
		gl dpvar = "delp"
		gl devar = "dele"  
		gl bvar = "budget"  
		gl stvar = "bstop"  
		gl sbvar = "bsbottom"
		gl ppvar = "pricep"  
		gl pevar = "pricee"  
		gl tvarp = "thetap"  
		gl tvare = "thetae"  
		gl mvar = "M"  
		gl np0var = "Np"  

		su Np 
		local p=r(mean)
		
		su Ne
		local e=r(mean)

		 mat b0 = ( `p' )
		 mat colnames b0 = $nlist
		 
		 amoeba find b0 yout bols . . 1E-10 
		 
*****************************************************************		

		 
		 
*****************************************************************		
** Calculate quantities and CS based on optimal holdings

		 
		gen mNp = bols[1,1]
		gen mNe = (budget - pricep*mNp)/pricee 
		
		su Np Ne mNp mNe 
		gen factor = 20
			gen double  Dgp = (mNp*exp(delp/(1-bsbottom)))^((1-bsbottom)/(1-bstop)) 
			gen double Dge  = (mNe*exp(dele/(1-bsbottom)))^((1-bsbottom)/(1-bstop))
			gen double sDg = Dgp + Dge 
			gen double mQp  = M*Dgp*(sDg^(-bstop))/(1+sDg^(1-bstop))
			gen double mQe  = M*Dge*(sDg^(-bstop))/(1+sDg^(1-bstop))
			gen double CS  = M*log(1+sDg^(1-bstop))
		
*****************************************************************		
	

*****************************************************************		
** save the variables we care about

		keep Np Ne mNp mNe lno factor Qe Qp  mQe mQp CS lno M dele delp thetae thetap ino
		
		gen raise ="Pp"
		
		sleep 500
		append using "$datapath/intermediate/collect_pdouble_income.dta" 
		sleep 1000
		capture save   "$datapath/intermediate/collect_pdouble_income.dta", replace  

	restore 
		
}

	
